home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / ode-initval / rk2.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  4.5 KB  |  205 lines
  1. /* ode-initval/rk2.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Runge-Kutta 2(3), Euler-Cauchy */
  21.  
  22. /* Author:  G. Jungman
  23.  */
  24. #include <config.h>
  25. #include <stdlib.h>
  26. #include <string.h>
  27. #include <gsl/gsl_errno.h>
  28. #include <gsl/gsl_odeiv.h>
  29.  
  30. #include "odeiv_util.h"
  31.  
  32. typedef struct
  33. {
  34.   double *k1;
  35.   double *k2;
  36.   double *k3;
  37.   double *ytmp;
  38. }
  39. rk2_state_t;
  40.  
  41. static void *
  42. rk2_alloc (size_t dim)
  43. {
  44.   rk2_state_t *state = (rk2_state_t *) malloc (sizeof (rk2_state_t));
  45.  
  46.   if (state == 0)
  47.     {
  48.       GSL_ERROR_NULL ("failed to allocate space for rk2_state", GSL_ENOMEM);
  49.     }
  50.  
  51.   state->k1 = (double *) malloc (dim * sizeof (double));
  52.  
  53.   if (state->k1 == 0)
  54.     {
  55.       free (state);
  56.       GSL_ERROR_NULL ("failed to allocate space for k1", GSL_ENOMEM);
  57.     }
  58.  
  59.   state->k2 = (double *) malloc (dim * sizeof (double));
  60.  
  61.   if (state->k2 == 0)
  62.     {
  63.       free (state->k1);
  64.       free (state);
  65.       GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
  66.     }
  67.  
  68.   state->k3 = (double *) malloc (dim * sizeof (double));
  69.  
  70.   if (state->k3 == 0)
  71.     {
  72.       free (state->k2);
  73.       free (state->k1);
  74.       free (state);
  75.       GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
  76.     }
  77.  
  78.   state->ytmp = (double *) malloc (dim * sizeof (double));
  79.  
  80.   if (state->ytmp == 0)
  81.     {
  82.       free (state->k3);
  83.       free (state->k2);
  84.       free (state->k1);
  85.       free (state);
  86.       GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
  87.     }
  88.  
  89.   return state;
  90. }
  91.  
  92.  
  93. static int
  94. rk2_apply (void *vstate,
  95.        size_t dim,
  96.        double t,
  97.        double h,
  98.        double y[],
  99.        double yerr[],
  100.        const double dydt_in[],
  101.        double dydt_out[], 
  102.            const gsl_odeiv_system * sys)
  103. {
  104.   rk2_state_t *state = (rk2_state_t *) vstate;
  105.  
  106.   size_t i;
  107.   int status = 0;
  108.  
  109.   double *const k1 = state->k1;
  110.   double *const k2 = state->k2;
  111.   double *const k3 = state->k3;
  112.   double *const ytmp = state->ytmp;
  113.  
  114.   /* k1 step */
  115.  
  116.   if (dydt_in != NULL)
  117.     {
  118.       DBL_MEMCPY (k1, dydt_in, dim);
  119.     }
  120.   else
  121.     {
  122.       int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1);
  123.       GSL_STATUS_UPDATE (&status, s);
  124.     }
  125.  
  126.   for (i = 0; i < dim; i++)
  127.     {
  128.       ytmp[i] = y[i] + 0.5 * h * k1[i];
  129.     }
  130.  
  131.   /* k2 step */
  132.   {
  133.     int s = GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h, ytmp, k2);
  134.     GSL_STATUS_UPDATE (&status, s);
  135.   }
  136.  
  137.   for (i = 0; i < dim; i++)
  138.     {
  139.       ytmp[i] = y[i] + h * (-k1[i] + 2.0 * k2[i]);
  140.     }
  141.  
  142.   /* k3 step */
  143.  
  144.   {
  145.     int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp, k3);
  146.     GSL_STATUS_UPDATE (&status, s);
  147.   }
  148.  
  149.   /* final sum and error estimate */
  150.   for (i = 0; i < dim; i++)
  151.     {
  152.       const double ksum3 = (k1[i] + 4.0 * k2[i] + k3[i]) / 6.0;
  153.       y[i] += h * ksum3;
  154.       yerr[i] = h * (k2[i] - ksum3);
  155.       if (dydt_out)
  156.     dydt_out[i] = ksum3;
  157.     }
  158.  
  159.   return status;
  160. }
  161.  
  162. static int
  163. rk2_reset (void *vstate, size_t dim)
  164. {
  165.   rk2_state_t *state = (rk2_state_t *) vstate;
  166.  
  167.   DBL_ZERO_MEMSET (state->k1, dim);
  168.   DBL_ZERO_MEMSET (state->k2, dim);
  169.   DBL_ZERO_MEMSET (state->k3, dim);
  170.   DBL_ZERO_MEMSET (state->ytmp, dim);
  171.  
  172.   return GSL_SUCCESS;
  173. }
  174.  
  175. static unsigned int
  176. rk2_order (void *vstate)
  177. {
  178.   rk2_state_t *state = (rk2_state_t *) vstate;
  179.   state = 0; /* prevent warnings about unused parameters */
  180.   return 2;
  181. }
  182.  
  183. static void
  184. rk2_free (void *vstate)
  185. {
  186.   rk2_state_t *state = (rk2_state_t *) vstate;
  187.   free (state->k1);
  188.   free (state->k2);
  189.   free (state->k3);
  190.   free (state->ytmp);
  191.   free (state);
  192. }
  193.  
  194. static const gsl_odeiv_step_type rk2_type = { "rk2",    /* name */
  195.   1,                /* can use dydt_in */
  196.   0,                /* gives exact dydt_out */
  197.   &rk2_alloc,
  198.   &rk2_apply,
  199.   &rk2_reset,
  200.   &rk2_order,
  201.   &rk2_free
  202. };
  203.  
  204. const gsl_odeiv_step_type *gsl_odeiv_step_rk2 = &rk2_type;
  205.